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Abstract 

Time scales have been constructed in different ways to meet the many demands placed upon 
them for time accuracy , frequency accuracy , long-term stability and robustness . Usually , no single 
time scale is optimum for all purposes. In the context of the impending availability of high-accuracy 
intermittently -operated cesium fountains , we reconsider the question of evaluating the accuracy of 
time scales which use an algorithm to span interruptions of the primary standard . 

We consider a broad class of calibration algorithms that can be evaluated and compared 
quantitatively for their accuracy in the presence of frequency drift and a full noise model (a 
mixture of white PM, flicker PM, white FM, flicker FM and random walk FM noise). We present 
the analytic techniques for computing the standard uncertainty for the full noise model and this 
class of calibration algorithms. The simplest algorithm is evaluated to find the average-frequency 
uncertainty arising from the noise of the cesium fountain’s local oscillator and from the noise of a 
hydrogen maser transfer-standard. This algorithm and known noise sources are shown to permit 
interlaboratory frequency transfer with a standard uncertainty of less than 10“ 15 for periods of 
30-100 days . 


Introduction: The Need for Evaluating Algorithm Accuracy 

For the near future, new primary (cesium fountain) frequency standards [1] [2] are likely to oper- 
ate intermittently, rather than to operate continuously as primary clocks. Other new frequency 
standards of high potential accuracy, such as single-ion optical frequency standards coupled to a 
divider chain [3], are also likely to operate intermittently, at least initially. More reliable secondary 
standards of high stability (such as hydrogen masers) will be needed to span the gaps between 
periods of operation of the primary standard. To evaluate the accuracy of time scales that are to 
be calibrated with these new standards, one must address the question of how the random noise of 
the primary standard, of its local oscillator, and of the secondary standards all combine to influ- 
ence the accuracy of the time scale or its average frequency. We want to examine how these noise 
sources affect the results of different interpolation and extrapolation algorithms, and to predict the 
accuracy that could be delivered, in the presence of mixed types of noise, to a local time scale or 
to TAI. Our main interest is in the frequency accuracy of the secondary time scale after calibration 


249 



using some algorithm, but many of the same ideas are also applicable to time accuracy around an 
interval. 

This is not a wholly new question. Intermittent operation of primary cesium frequency standards 
was universal two decades ago. Continuous operation of these standards as primary clocks, first at 
NRC [4] and then at PTB, waited until time laboratories had sufficiently evaluated their primary 
cesium frequency standards and had improved their mean time between failures. Standard tech- 
niques for analyzing frequency standards’ stability and characterizing mixed types of noise were 
adopted and refined [5] [6] [7], A body of very useful guidance was developed [8] for extrapolation 
in the presence of different (but unmixed) types of noise. 

The main new element to be addressed is the quantitative estimation of accuracy for the mixed 
types of noise which have to be faced for our problem, and which has not been needed for previous 
standards where the dominant noise has usually been flicker frequency noise for primary and sec- 
ondary standards. Simulations [9] could give the required guidance, but fully analytic techniques 
are preferable, and are developed here for a widely used class of algorithms. 


Metrics 

In choosing and in judging accuracy of the “optimum” algorithm for a purpose, a metric should 
be used for ordering possible algorithms and for guidance of minimal ambiguity. A priori, there are 
many possible metrics. 

The class of metrics of interest to us quantifies the difference between any two functions of time, 
A(t) and B(t ), sampled at a set of discrete times {£j},i = 0..N during the time interval [£o»*iv]- A 
metric expresses as a real number the difference between two vectors A and B in this iV-dimensional 
vector space: l|A-B||, and permits the unambiguous ordering of the quality of a fit from “good” 
to “bad”. Any metric must meet the requirements that (1) ||A|| > 0, (2) ||aA|| = |a|||A||, and 
(3) IIA + BH < || A || + ||B||. A very useful subclass of metrics is the class of “Holder norms”, 
or L P -norms (p > 1): ||A — B|| = E^Lo | (A* — Bi ) The weights Wi are positive 

definite. If the difference between A{ and Bi is a random variable z x distributed around Z*, and 
is described by a probability distribution then the minimizing the corresponding 

L p - norm will give the maximum likelihood fit of A to B. When fitting an approximation to a 
mathematical function, the norm (limp_ +00 ) is usually used, as the min-max norm, to minimize the 
maximum error between the function and its approximation on the interval. The absolute- value 
norm (p = 1) is occasionally used as an uncritical way of fitting to give minimum fitting weight to 
erroneous outliers while formally retaining a metric. When fitting experimental data, where normal 
(gaussian) distributions are common, p=2 is generally appropriate. It is appropriate for describing 
our expected distributions, and we will concentrate on this type of metric. 

When measuring the quality of a fit to the measurements at the N times £i, the value of the metric, 
divided by the degrees of freedom ( N minus the number of fitting parameters), is often used. For 
a least-squares fit (p = 2) this measure of the quality of fit becomes the square root of the familiar 
reduced x 2 i and for unweighted least-squares fits ( Wi = 1) it is the even more familiar root-mean- 
square residual. The residual is formally a metric in an N-dimensional vector space. As they are 
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formally defined, weights can in principle be used to change the scale factor on each axis of the 
vector space, and even to project the metric into a vector subspace. Minimizing such a reduced 
metric rather than using the full maximum-likelihood weights can be advantageous: for example, 
the 10-day average frequency of a commercial cesium clock can be better determined from only the 
end points than by a least-squares fit to many points distributed across the 10-day interval [8], 

Thus the quality of the fit at the experimental points is not all that is required. Rather, since the 
fitted function is used for interpolation to other values of t between U and ti+ 1 , or for extrapolation 
to values of t outside [£o>£jv]> it is the accuracy at these points which can be more important. 
The accuracy of an algorithm is not uniform, but varies with t in a way which depends on the 
set of fitting points {U\i = 0..7V}, the fitting algorithm and the random (and deterministic) noise. 
Considering this type of problem from the perspective of the residuals seems to require the magic of 
rotations into a different vector space. Interestingly, exactly this task can be done for the L 2 norm 
and a rather broad class of fitting functions, although the metric projection picture is unhelpful in 
determining the fitting accuracy at an arbitrary time. 

.The accuracy can be determined for any system experimentally by repeated cycles of measurements, 
doing repeated fitting of one particular pattern of time samples and by statistical analysis of 
residuals determined at unfitted points. Another approach would be computer simulation of this 
process - if a sufficiently good description of the noise model is available; or it might be done 
analytically. We show that a rather broad class of noise models and fitting procedures can be 
treated analytically, to obtain an accuracy estimate for the interpolation or extrapolation of many 
commonly used algorithms. 


Modelling the Difference: Deterministic plus Random Noise 


To describe the time-dependence x(t) of the time difference between the primary frequency stan- 
dard’s time scale and the secondary time scale, we model it with x m (t) ) and explicitly include a 
random part xq ( t) as well as a deterministic part. The deterministic part allows for a time offset, 
a frequency difference, and a drift rate of the frequency of the secondary time scale with respect to 
the primary standard. 


x m (t) = a k + b k t + ^-c k t 2 +x 0 (t). (1) 

z 

The superscript k labels the uninterrupted intervals of operation of the primary standard. For each 
interval, a new value of a k is required, and other values of b k and c* may (or may not) be used. 
We will examine the accuracy of a class of fitting functions x p (t ), fit on the interval A;, linear in the 
fitting parameters {d/}, and including a broad class of basis functions g\{t ): 


x p(t) — d k + d k t + d k t 2 + ^2 dj[9i(t). ( 2 ) 

1=4 

In the random noise part, xo(t), we include the “full” noise model that is usually covered in discus- 
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sions of frequency standard stability [5]: a sum of five noise processes, each normally distributed 
about the mean (but with variances which depend on the time sampled in different ways) that have 
spectral densities of phase noise ( S x (t )) that are power laws which range from flat to increasingly 
divergent at low frequencies. Expressing the five terms in terms of the spectral density of the 
mean-square of the fluctuations in (or yo(t)) at a frequency /, S y (/), each noise term is de- 
scribed by an amplitude h a which is taken to be independent of any time translations (stationarity 
and random phase approximations). The sum includes (a = 2) white phase noise in x, (a = 1) 
flicker (1//) noise in x, (a = 0) white frequency noise and random walk phase noise, (a = — 1) 
flicker frequency noise and (a = — 2) random walk frequency noise. The spectral density of the 
mean-square fluctuations in xo(£) is S x (t), and for this noise model 

S y (f) = EL -2 h a f a and S x (f) = EL- 2 * 1 ^- (3) 


A Useful Tool 

For many cases we might expect mean-square accuracy estimates to be calculable from the auto- 
correlation function (xo(£)xo(£ + r)). Although it is not easy to use, the autocorrelation function of 
frequency-standard noise has been rather neglected in our view. It is divergent for four of our five 
types of noise unless a low-frequency cutoff is applied, and even then can challenge the accuracy 
and dynamic range capacities of classical computing. Analytic expressions for this autocorrelation 
function are given in the Appendix for each type of noise, and modern arbitrary-precision com- 
puter languages should routinely be able to cope directly with the autocorrelation function. In our 
analysis of the uncertainty associated with any useful time algorithm we expect no divergences to 
infinity, and so the combinations of the autocorrelation functions must have their divergent parts 
cancel, with the algorithm itself acting as low-frequency cutoff. To simplify the numerical evaluar 
tion of combinations of the autocorrelation, it can be useful to find an analytic expression for the 
general two- interval covariance of the random noise model, that is the covariance of the time-scale 
departure over the time interval with the time-scale departure over the time interval [£3, £4]: 

<5 = ([xo(t 2 ) - xo(ii)] M4) - x Q (t 3 )]) = Su(yo(t') yo(t ")) dt " dt '- (4) 

This is just (t 2 — ti)(tn — £3) times (5[t ll t J ]&[t 1 ,t 4 ])j the general covariance of the average frequency, 
which is a generalization of the two-sample variance of the average frequency. The generalization 
includes the possibility of an overlap of the intervals (as well as the possibility of a “dead time” 
between the intervals), and incorporates the possibility of considering the frequency average over 
two time intervals of different duration. Just as for the two-sample variance of y, and for the 
autocorrelation function of x(£), the covariance separates into the five terms of the noise model. 

Analytic forms for the five terms of the autocorrelation function of x(£) and for the five terms of 
the general cross-correlation of y are presented in as Equations 18 to 25 in the Appendix, derived 
with only the usual assumptions about high and low frequency limits to the noise bandwidth. Some 
comments are made on practical methods for computing values using these forms. 
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Algorithms for Accuracy Evaluation 


For many, but not all, widely- used interpolation or extrapolation algorithms, it is possible to use 
the cross-correlation expressions and a knowledge of the noise to calculate the expected root-mean- 
square (i.e. p = 2) standard uncertainty [10] of the time or of the mean frequency extrapolated 
or interpolated by the algorithm. Based on the noise model (and a very large body of confirming 
experiment), the distributions expected for deviations from the fit are normal, and so the root- 
mean-square deviation calculated for the standard uncertainty is rigorously correct, and can be 
used in the conventional way for predicting confidence intervals from normal distributions [10]. 
The fitting process yields parameters for a parameterized functional form of x p (t), which may be of 
the general form (around the k ih interval of continuous operation) given in Equation 2. Nonlinear 
fitting parameters are not considered. 

The algorithm should satisfy two criteria. Firstly , the algorithm should be unbiased by the deter- 
ministic part of the noise model: any change in a deterministic parameter (a fc , 6* or c k in Equation 
1) should be taken up by the algorithm and not bias the final result. Note that it is the final 
deviation which is to be unbiased, and some apparently biased forms for x p (t) may still be used in 
ways that are unbiased. In addition to being patently desirable, this condition also removes any 
need for considering cross-correlations between the deterministic and random noise parts of x m (t). 
Secondly , the coefficients d\ must depend linearly on sums over differences in x(t). This includes 
fitting functions that are constrained to go through one point: least squares fits of polynomials of 
general order, and other functions with linear coefficients. It includes constrained weighted least- 
squares fits, as long as the weights do not themselves depend on the values of x k (ti) or the variances 
on the k th fitting interval. With this condition we ensure that we do not have to calculate any 
higher-order correlations than the general covariances evaluated in the Appendix. 

What algorithms does this exclude? The first condition would not seem to exclude any serious 
contenders for fitting methods and fitting functions: if one is taking care to evaluate accuracy, then 
presumably one also values an unbiased fitting form. The second criterion admits many algorithms 
easily. However, to analyze rigorously the accuracy of a Kalman filter, where fitting weights depend 
on past variances, appears to require a study of higher-order autocorrelations of the random noise, 
at least to the level of identifying the magnitude of these corrections. We therefore exclude this 
important class of algorithms from our present considerations. 

Many extrapolation and interpolation algorithms that are useful for time-scale purposes are very 
simple: for example, constrained to go through the last experimental point [8], or constrained to 
go through both the first and last points on an interval. However, it is instructive to consider first 
the most general least-squares fitting case for N points on a single calibrating interval. 

Weighted Fits of General Functions with Linear Coefficients 

The least-squares fit of the n parameters di of a function of the form of Equation 2, to N + 1 
points x(U), each point being taken with a weight u;*, is found by taking the partial derivative of 
the weighted L 2 norm with respect to its n coefficients. The resulting n simultaneous equations 
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in n unknowns are of the form G 3 = s', where G is an n x n matrix with elements G qr = 
Y^iLo w i9q(U)9r(U)) and & is an n-dimensional vector with elements s' r = w i x m(U)9r(U)- Here 
x m (ti) is used to model x(ti ), exactly as one might do in a numerical simulation. Simple inspection 
shows that if g\(t) = 1,02(0 = t and g^(t) — t 2 , and s r = w i x o(U)9r(U)i then Gd = s , 
where d is the vector of least-squares coefficients for the n functions, and d\ = d[ — a, d<i = d f 2 — b , 
cfe = d^ — c/2 and d*. = dj. for r > 4. Thus for any set of weights, the fit exactly absorbs the 
deterministic part of the model function x m (t) and we have only to deal with the function xo(t). 
We can replace any difference between x m (t) and x p (t) (fit to the x m (t*)) with the difference between 
xo (t) and the x p (t) (fit to the xo(t»)). The fitting coefficients vector d = G _1 s*. The form of this 
equation is instantly quite revealing, for it shows that the least squares coefficients dr only involve 
simple sums over the xo(^i)’s. 

We will study all the effects of algorithm at t , reacting to the noise through the fitting procedure, 
by comparing xo (t) to the function fitted to the random part of the noise: d • g(t ), where g(t) is 
the n-dimensional vector of the basis functions evaluated at t. This fitted function can in turn be 
transformed into a simple weighted sum over the N + 1 of the xo(£i) 5 s: d* g(t) = A(O x o(£i), 
and Di(t) = (U)g q (t). 

The fitting algorithm and the noise model contribute a standard uncertainty [10] u x in the least- 
squares fitted time which is ([x m (t) — x p (t)] 2 ). This is the appropriate metric for judging the quality 
of the fit at t (relative to the set of N + 1 specific fitting times {£;} ). 


([x m (0 -Xp(£)] 2 ^ = ( [x m (*) -(?•(?(£)] = ([x 0 (£) -<T- <?(£)] 

when t is labelled as and D-\ = 1 in the last equation. The sum over the (N + 2) 2 autocorrela- 
tions simplifies since the autocorrelation depends on | t{ — tj |. When the data are equally spaced, 
“only” N + 2 different autocorrelations of xq ( t) must be evaluated for a general value of t. The 
autocorrelations for our noise model can be evaluated using the expressions in the Appendix for 
1 £(t), if calculations can be done with sufficient numerical precision. 


Constrained Least-Squares Fits 


Fits constrained to go through one or more points can be considered as special cases in the above 
analysis through appropriate choices of weights. However, there is an interesting computational 
simplification which warrants explicit derivation: we can replace computations of the autocorre- 
lation function of x(t) with the simpler to compute S of Equation 13. Consider using a weighted 
least-squares fit constrained to go through one of the x m (^)’s, for i = i c in general - although this 
might often be either end point: i = 0 or i = N since we consider {£;} as being ordered so that 
tj < tj+ 1 . For a constrained fit, x m (ti c ) = d f c • g(U c ), determined by substituting g(U c ) for the first 
row of G to obtain G c , and by substituting x m (ti c ) for the first element of s' to obtain ^ c . The 
weight Wi c is assigned the value 0. The constrained fit is also unbiased by the deterministic part 
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of the noise model, satisfying the condition so that again we only have to consider the behaviour 
of xo(t) when evaluating the accuracy of the constrained fit. The fitted function is d c ■ g(t ), where 
d c = and s c is just with xq ( t) substituting for x m (t). For this constrained fit we can use 

the equality of x m and x v at or the equality x$ (ti c ) = 0: 


-2Tp(0] 2 ) = ^[{®m(0- a! «( t J}-{ < £-(^(0-fl( t i e )}] y = 

[{x 0 (t) - zota)} 

■ {Eg=l {9q(t) — 9q(U c )} {(G^ 1 )^ Xo(ti c ) + Er= 2 (G c l )qrJliLo w i9r(ti) ^o(^i)}}] ^ = 

[{xo(0 -X 0 (t,J} - {Zilo^cM^}] 2 ) = 

[{*o(0 - Z0 (tic)} - {Zlo(D c )r {E }=1 (*0 (tj) - *o(t,--l)) - E }=1 (xo (tj) - x 0 (t,-l))}}]‘ 
= ([Wt) - X 0 (t lc )} - {Ef=l(4); {*o(f;) - ^o(t;-l)}}] : 

= ([{EjLo(4);(*0(t;)-®o(ti-l))}‘ 


( 6 ) 


where (D c )j is just the reordered sum defined above for j — 1 to iV, with ( D c )q defined as 1. 

Thus the constrained fit’s standard uncertainty in time can be calculated more easily at a general 
time t, using only our expressions in the Appendix for I(r) rather than for TZ x (r) to evaluate the 
standard uncertainty in time, which can done conveniently with 64-bit floating point calculations. 
Note that this is also the mean square of the time interval error over the time interval [t lc , t]. 


Standard Uncertainty in Average Frequency 


In a similar way we can calculate our model’s estimate for the unconstrained least- squares fit’s 
standard uncertainty [10] of the average frequency over an interval [t,t -f r], u v (t,r), We can 
calculate u 2 (t,r) = + r) - x m (t)) - (x p (t + r) - x p (t))} 2 } /r 2 in terms of expressions using 

only the function T(T). For convenience, we will define t_i = t and tw+\ = t + r, neither restricting 
the value of t to be less than to nor restricting the value of t + r to be greater than t;v. 

u y( r ) r2 = +T) - X m (t)) - (X p (t + T) - X p (t))] 2 } = 

l [{^o(£ + r) - Xo(0} - {d- (§(t + T ) - 9(0)}] ^ = 

([Ef=VN(t;)-^-i)} + Ef=i^ {^o(ty) -a:oOy-i)}] 2 ) = (?) 

( [E^o 1 D } {x 0 (t ; ) -x 0 (t ; _i)}] 2 ^ 

where Dj = T.tj v>f EJ=i E”=i(G -%9r(ti)0h(t + r)~ 9,(0) and Dj = [Dj + 1] for j = 0 to N, 
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and where both Do and Dn+\ are equal 1. Again the uncertainty can be calculated conveniently 
with 64-bit floating point calculations from the expressions in the Appendix for I(r) rather than 
the more awkward 72 x (r). 

Multiple Calibration Runs 

The time interval error that accumulates on an interval spanning m calibration runs is just the sum 
of time interval errors (the difference between the evolution of x and the evolution of the calibrated 
extrapolation x p ) on the interval around the k th calibration run. The calibrated time scale 

is continuous, so that if the k th fit yields to the next fit (fc + 1) at a time t r * = then 3* ■ g k (t r k ) = 

d fc+1 -S* +1 (t| *+ 1 ). Thus, across the m runs, the time interval error of the algorithm reacting to our 
noise model is 

m m 

£ J5* = £ M‘r‘) - *°(V)] - kM - *}(«!•) ■ (8) 

k= l fc=l 

To show the general form, consider the set of {t^}, with the index r) running in turn over the start 
time, the fitting times and the stop time for each interval, from 1 to M - the grand total of points 
(with N k measurable subintervals and two extrapolated subintervals in the k th interval). To analyze 
a fit on the k th interval that extracts information from other intervals, we consider the k th fit to 
be from t L k to t R k , with zero weight at the unmeasured times t v where r) = V or r- 7 , j = 1 . . . m. 
Equation 8 can then be combined over any overlap of the fits into a global weighted sum over the 
differences [xq ( t v ) - xo(t v -i)]> weighted by T> v which sums over the fits which have used the rf h 
time interval. The mean square time uncertainty in the time scale algorithm over the total time 
interval T = Y^f= iftv ~ *i?-l] 13 the mean square of the sum over Ei . Since for any algorithm 

of the class 


j ( M \ 2 \ /MM \ MM 

( =EE (wg - x <>( vi)i fade) - ®o(t f -i)]> . (9) 

\v?=i / ! c=i 

The full evaluation of all terms of this M x M cross-correlation matrix, £ (where £ij = EiEj ) 
is possible for any particular set of calibration runs. With M 2 terms to evaluate, an efficient 
method for computing S is highly desirable. The method outlined in the Appendix will generally 
suffice. Fortunately, there can be very significant simplifications: £ is symmetric, the main diag- 
onal contributes most to the sum, the block-diagonal terms from the individual fits are the next 
most important parts (together these would contribute the quadrature sum of the standard time 
uncertainty u x contributed across each individual interval, but neglecting interval-to-interval corre- 
lations), and generally the matrix elements far from the diagonal will not contribute significantly. 
Any regularity in the fitting times will also reduce the computational burden. 
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Examples 


We will apply the above methods to one of the simplest algorithms that might be used: a linear 
fit constrained through two points, as shown in Figure 1. Here the mean square of the standard 
uncertainty in time is [11] 


t £=2 


(1 + —)T{t) + — (1 + —)I(t a ) 


ta 


I(t a + T)\ 


( 10 ) 


This expression can be separated into the standard uncertainty u x for each noise type, as illustrated 
in Figure 2. In this figure, u x for each noise type has been normalized to equal 1 at the midpoint 
of the calibration interval t a . The upper frequency cutoff u> is 100/i a , and the lower cutoff e = 
10~ 5 radians/second. To use a figure like this for a mixed noise model, recall that the standard 
uncertainties of the different noise types must be added in quadrature, with appropriate weights. 

The above example shows how an rms residual at a time offset r can be calculated, and that even 
in the simplest case it varies with position in different ways for the different types of noise. Using 
X(r) it is calculable in a perfectly straightforward manner. 

Figure 3 shows the time interval error that develops from extrapolation to both earlier and later 
times than the calibration interval: from £/ earlier than the first point and to t r later than the 
second fitted point. The standard uncertainty in frequency, u y (r), on this extended interval of r is 
[14] 


U 2 y = ^|l( T )+(^) 2 I(t a ) + ^I(l r )-^ 1 + i a ) + ^)-^a + ^)} (11) 

As an interesting application of this simple algorithm, we can show how the local oscillator limit 
[12] might be circumvented for a pulsed cesium fountain. Consider the evaluation of different types 
of local oscillator for a pulsed cesium fountain that employs this algorithm and hyperfine phase 
differences [11] to span 0.01 s dead times between 0.99 s Ramsey times. The atomic polarization 
is measured for each pulse of atoms, and attributed (after calibration) to the discrepancy of the 
microwave phase compared to the primary hyperfine phase of the ensemble of cesium atoms. Thus 
the frequency of the local oscillator is known across the 0.99 s interval, with an uncertainty that may 
be limited by (among other things) atom shot noise in the ensemble. The average frequency from 
the preceeding and following Ramsey time is used to estimate the frequency of the local oscillator 
across the 0.01 s dead time between Ramsey pulses. If there is but a single ensemble “in flight” 
through the cesium fountain at any given time, then the simplest algorithm is the only one worth 
considering. 

To analyze this system, we choose a symmetric interval centred on t a) the active time, and 1 1 = 
t r = t d j 2. Using the calibrated frequency from a single active interval for the two adjacent dead 
time half-intervals is equivalent to using the average frequency from two adjacent active intervals 
across any dead time. In this way we might hope to minimize the cross-correlations that need to be 
evaluated between neighbouring active times. We can use Equation 11 for the first estimate, but we 
should verify the size of the correction from neighbouring intervals. With a cycle time t a + t d = T, 
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and N equal-length cycles, Equation 9 becomes 


£ j+lij = 2{IT + T)+ 2{IT - T) - 2(1 + [£] 2 )Z(iT) + [£] 2 {I(IT + t a ) + 2{IT - t a )} 

-2J {I {IT + \[T + t a ]) +2 {IT - k[T + t a }) -2{IT + \{T - g) —2 {IT — ^[T - g)} . 

Fortunately, for l > 1, the off-diagonal corrections are negligible. Depending on a, the adjoining 
term, l = 1 , gives corrections to the diagonal terms of either sign and of up to 50% in magnitude. 
Thus correlations are significant in this case even though we tried apply our algorithm in a way 
to minimize the correlations. Figure 4 shows the standard uncertainty in frequency due to the 
local oscillator noise for three local oscillators. The curvature reflects some of the effects of cor- 
relations. The Allan deviation predicted [12] for the same three local oscillators in a conventional 
frequency servo is also shown (with a modulation frequency of 0.5 Hz). Clearly this hyperfine phase 
measurement plus post-processing algorithm is an interesting possibility as a replacement of the 
conventional servo, at least in some applications. It is crucial to have full confidence in the analysis 
of the complete effects of the interaction of the noise and the algorithm: our analysis is complete 
within the constraints of the noise model. 

In the initial period of operation of a cesium fountain frequency standard, after having been eval- 
uated as a standard, there might be a hour per day devoted to calibrating a time laboratory’s 
secondary standards. In the case of NRC, this would mean the use of our new hydrogen masers 
[14], [15]. A crucial question in designing a cesium- fountain frequency-standard project is to choose 
the desired level of accuracy which might be used, either transferred to another laboratory or used 
to compare the frequency of two configurations of the cesium fountain which cannot coexist. The 
answer to this question can be obtained by our method for secondary frequency standards with well 
understood noise. As always, common mode noise between similar frequency standards is difficult 
to rule out - but the intent of these considerations is to establish a goal for an initial working 
standard that is not overbuilt, considering available frequency transfer characteristics. We present 
the analysis of a possible frequency-transfer budget at NRC. 

The main question is the frequency transfer capabilities of our hydrogen masers. The Allan devi- 
ation of our two new masers has been measured with respect to each other. Attributing the noise 
equally to the two masers, we can calculate the best case for frequency transfer from one hour 
out to a day or so [14]. We again use the symmetric linear extrapolation from the end points. In 
this example, additional information will be available during the interval, but the “best case” for 
our noise types comes from using the end points [8], By using the symmetric extrapolation, the 
frequency transfer will not be biased by any constant drift of the maser frequency [13], 

Figure 5 shows an estimate for the standard uncertainty in frequency due to a hypothetical cesium 
fountain [11], [15], and the modelled Allan deviation one of the new NRC hydrogen masers [14]. 
The random-walk FM rise at large times is a somewhat pessimistic (or realistic ?) guess - the 
masers have not operated for a long enough time to properly evaluate their long-term stability 
<t v (t > 30 days). Also shown is the calculated standard uncertainty of the average frequency of the 
fountain-plus-maser for the extrapolated frequency on an interval r for a calibration run of 3000 s. 

Using Equation 12, we can again evaluate cross-correlations for a regular series of runs (daily in 
this example). However, in this example there are no corrections larger than 1% to the N ~ 1 / 2 trend 
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line. At each end of r we include a two-way time transfer uncertainty of 2 ns (1 - a per transfer) 
as a realistic estimate of the state-of-the art (and including 0.2 ns as a worthwhile objective). The 
daily recalibration trend line gives the noise limit to the standard uncertainty for cesium fountain 
frequency transfer to another laboratory. It is encouraging in that the noise limit is well below 
10" 15 for an ambitious but realistic cesium fountain calibration schedule for the hydrogen masers. 
The estimated standard uncertainties for the as- transferred average frequency may facilitate the 
acceptance of the new standards, for example for the judging or steering of TAI. The predicted 
level of residual fluctuations will need to be experimentally verified to be fully confident that some 
significant common-mode noise source (random or deterministic) has not been missed, but even 
here having a baseline prediction will be very helpful in planning the evaluation level which is 
appropriate for any given frequency-transfer program. 


Conclusions 

We have developed, calculated and applied some useful metrics for judging the accuracy of algo- 
rithms extrapolating time and average frequency in the presence of noise that can be represented 
by a rather general noise model which includes all common types of power-law random noise as well 
as deterministic noise. For many algorithms the (rms) standard uncertainty in time, u x (t), and the 
(rms) standard uncertainty in average frequency across r, Uy(t,r), both can be calculated in terms 
of the noise model power law coefficients. This significantly enhances the attractiveness of standard 
uncertainties in time and frequency metrology where techniques for measuring the coefficients are 
widely used. The metrics can be used for guidance in the choice of algorithms and procedures (how 
often to recalibrate, and for how long), but a larger role can be played by these two “metrics”, for 
rigorously judging the noise floor of different hardware and potential hardware used in novel ways. 

The calculated standard uncertainty in time, u x (t ), might also play a very useful role in the calcu- 
lation of a reduced x 2 for a set of experimental data for which one wishes to judge the adequacy of 
the fit and noise model. Conventionally, x 2 = 7 ^ £.=1 [*(*0 - z P (t>)] 2 /{2[u x (f,)] 2 }, where N - n 
is used for the degrees of freedom: the number of (independent) data points minus the number 
of (independent) fitting parameters. Since we would be able to compute the cross-correlations 
between data points, it should also be possible to determine a better estimate of the degrees of 
freedom for the fit. 

The procedures outlined here could be automated. Standard uncertainties could be used routinely 
in time and frequency metrology, for many common fits to a set of data points, in the presence of 
a general noise model. In all fields of metrology, where the use of standard uncertainties is now 
recommended [10], the very valuable techniques (such as the use of the modified Allan deviation) 
developed for analysing power law noise of frequency standards could be applied in other fields to 
give rigorous results for the standard way of reporting uncertainties. 
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Appendix: Autocorrelation Functions and General Interval Co- 
variance for Power- law Noise Spectra 


Consider two time intervals that may (or may not) overlap, and are not necessarily of equal duration. 
The general interval covariance of the random part of the time scale difference accumulated on the 
time interval [ti, *2] with that accumulated on [£3, £4] is 

5 = ([ao(*2)-*o(ii)][a:o(t4)-*o(t3)]) = <^o(«2)a:o(i4)> + <^o(ti)a:o(«3)) 

- ~ (xo(h)x 0 (t 4 )) . ^ ' 

where {xo(t)xo(t — r)) is the autocorrelation function: TZ x (r) = / 0 °°S X (/) cos27r frdf. We use 
the usual [5] upper-frequency limit f u = and low-frequency limit fi = e/( 2 ir). The sharp 

upper cutoff is an artifice, although one which could be implemented with a digital filter applied 
to the output of a phase comparator. The lower frequency needs to be low enough so that it 
does not perturb the low-frequency rolloff supplied by the fitting function. A too-low value for e 
will exacerbate the numerical precision problems in calculating and using 7?. x (£), which diverges as 

e — ► 0. For the usual noise model of S x (f) = Yli=- 2 > we ^ ave 

2 h r0 _a ) f UT 2 

n x( T ) = Y 7o~w +T / 2 cos u du — Y R a(r). (14) 

a =-2 JtT a =—2 

It is possible to express S as a sum over the functions /^(r^s, or as a sum over functions of similar 
form / a (r)’s that are less divergent as e — > 0: 


— [JZ x (t4 — £]) + ~ £2) — ^(£4 ~ £2) — ^ x (£3 “ £1)] 

T (£4 — £1) + T (£3 — £2) — ^(£4 — £2) — I (£3 — £j), 


(15) 


where I(r) is 


2 h r(l -a ) f UJT ^ 

J ( T ) = Y /o-Aa +r / -COSu}<fci= Y / «( T )- ( 16 ) 

a — -2 JtT a =- 2 

With the help of mathematical tables or a symbolic algebra program such as Maple V, the integrals 
for R a (r) and I a (r) can be done analytically for our values of a. 



are the integrals for white phase noise (a = 2). I2 is proportional to the high frequency cutoff 
as is the Allan deviation. 
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(19) 


fliCO = (0y * Ci ^ ~ Ci ( eT ^ 

I] ( r ) = t^t 2 t'y + ln ( wr ) _ c i ( wr )} ( 2 °) 

are the integrals for flicker phase noise (a — 1). Ci(x) is the Cosine Integral function, which can 
be easily approximated by using numerical approximation as in Abramovitz [16] for arguments 
greater than 1, or by direct numerical integration of the expression between brackets (see end of 
this Appendix). The time interval r tempers the logarithmic dependence on the frequency cutoff. 

Ro(r ) = --±t{ i—L + Si(ur) — - 5t(er) \ (21) 

27 r 1 ur er J 

I 

27T l U)T ) 

are the integrals for white frequency noise (a = 0). Si(x) is the Sine Integral function, which can be 
treated in the same way as the Cosine Integral for arguments greater than 1, or by direct numerical 
integration of the definition of Si(x) for small arguments. This term depends linearly only on the 
time interval r, and not on value of the high frequency cutoff. 


= fc-iy { 


r 2 f ( — cos (ujt) sin(cjr) 


— Ci(ujr) > — 


— cos (er) sin (er) 
_ — T 


-Ci(er) | (23) 


/_i(r) = h-i— 1.5 - 7 - ln 


( er ) “ | 


1 — cos (cjt) sin (cjt) 

/ To + 


Ci(cjr)| 


are the integrals for flicker frequency noise (a = —1), where 7 is 0.57721..., Euler’s constant. 
The logarithmic term looks as if it will diverge to infinity for a low frequency cutoff as e — ► 0, 
but in combinations like in S normally the combination of terms using /_i(r) is such that the 
terms multiplying 1.5 — 7 will cancel out, as will the logarithmic divergence with e, leaving a term 
depending only on the square of the time interval and a geometric structure factor that depends 
on the logarithm of ratios of the time intervals concerned: for the Allan deviation this is ln(2) as 
expected, and for 5 it is ln 


= 7ih 


Sin(€T) . COS [tTj 

T7F + 


+ ™ + Si(u>T) 
+ s Her) - 2 


r 3 [ 3 \ l-cos(wr) sin (wr) cos(wt) 

/_ 2 (r) = rt_ 2 y --{2 , .3 TutV ujr ^ 
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are the integrals for random walk frequency noise (a = — 2 ) for a low frequency cutoff e —► 0 . 
As expected, this expression will diverge to infinity in the limit e — ► 0 , but the combination of 
terms will cancel that divergence in many of our accuracy prediction problems, since our fitting 
procedures will act as an effective low-frequency rolloff for the residuals of the fit. In practice, 
the low frequency cutoff will depend on the observation (active) time and is always some value 
greater than zero. In our analysis of uncertainty, we should use the same value of e that was used 
in determining /i _2 from measured Allan variances of the standard in question. 

Using standard expansions for the Sine and Cosine Integrals [16], it is easy to compute with 
enough accuracy the values needed for arguments greater than 1. For arguments smaller than 
one, numerical integration can be done easily for Si(x). For Ci(x), the following transformation is 
helpful, since the second integral is easy to do numerically for small arguments: 

Ci(x) = 7 + In x - /* l -=f^du = 7 + In * - f 0 x ^ du . (27) 

In debugging any code written to evaluate 7Z x (r) or J(r), it is worthwhile noting that the two- 
sample variance or Allan variance is 

(T 2 v {t) = = 4X(r)-J ( 2Tj ' where (28) 


n x { 0 ) = 


hi 

(2tt) : 


-(w - e) + 


h i 

(2tt) : 


In 


7 + S<‘" 


■ U) 


-1 


) + -^7-( e 


-2 


• U 


- 2 > + 


Ld 


which should reflect the traditional expressions [5] when the typographical errors have been cor- 
rected (their equations 101 through 105 have had the term 7 + ln(27r// l r), where 7 is Euler’s 
constant, incorrectly typeset as 2 + ln(27r f^r) ). Some approximations used to obtain their tradi- 
tional expressions have not been made. The expressions given above for /«(r) are correct for all a 
even for cj c t < 1 , so that differences are to be expected if calculations are made for r near to or 
less than l/(27r/ c ), where f c is the upper frequency cutoff or bandwidth of the phase measurement 
system used to measure x(t): the expressions above give the correct results. 
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Figure 1. A simple algorithm: linear extrapolation or interpolation from a line constrained to lie through two 
points on x(t), separated by a time t a . Its standard uncertainty in time, u x (r), can vary with extrapolation time T, as 
shown in Figure 2. 



Figure 2. The standard uncertainty in time, u x (r), developed by the algorithm of Figure 1 for an extrapolation time t. 
It is shown separately for the five types of noise, each normalized to 1 at the midpoint of the fitting interval ta. In 
this example, (Oh = 100/t a . 
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Figure 3. A simple algorithm: linear fit constrained to lie through two points on x(t). The variation in frequency 
from the calibration interval t a is illustrated. With ti = t r , this is the basis for Figures 4 and 5. 



io° io 1 io 2 T io 3 io 4 

Figure 4. Local oscillator contribution to the standard uncertainty of the average frequency for a pulsed-ensemble 
cesium fountain, with a cycle time of 1 s and a dead time of 0.01 s. The light lines show the classical stability limit 
of the three oscillators, and the heavy symbols show the pulsed-ensemble result using linear extrapolation in phase 
to bridge the dead time. For each type of servo, the top curve is for an Oscilloquartz 8600-3 (h-i =8xl0' 26 , 
hi = 8xl0' 27 ,h2 = 5.6xl0- 29 ), the middle curvefora Wenzel 500-03475 100MHz&5MHz(h.i =8xl0 26 , hi = 1.6X10' 30 , 
h2 = 1.3xl0' 34 ),and the bottom curve is for a JPL-type77K sapphire X-band frequency discriminator (h-i = lxlO" 27 ). 
Optimally used, their noise corresponds to a shot noise of 700, 8xl0 4 and 6xl0 6 atoms/s respectively. 
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t (seconds) 


Figure 5. Cesium fountain’s standard uncertainty of average frequency over an interval t,due to random noise. 
The maser noise model ’s Allan deviation is also shown as the heavy curve. A cesium fountain with a u y (t) = 10~ 14 T 1/2 
is operated for 3000 s , calibrates a maser, the standard uncertainty of the extrapolation from one calibration is 
indicated by the light curve that rises abruptly at t = 3000 s . The other curves show what can be done with current 
(2 ns) and a possible future two-way time transfer at the ends of the interval X . The dots show what can be done if a 
3000 s calibration run is performed every day. 
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